! **************************************************************************************
!  Rob Crain (2009), Centre for Astrophysics & Supercomputing, Swinburne University
! **************************************************************************************
! Code to identify all particles in P-Gadget3-BG snapshots that need to be
! resampled in a hi-res mask when generating intial conditions for resimulations.
! **************************************************************************************

program SelectParticles
  use log
  use myfiles
  use particles
  use params_and_constants

  implicit none
  integer*4 :: isnap

  ! Initialise log files and say hello
  call startup()

  ! As it says on the tin...
  call parse_paramfile()

  if(mode < 2) then 

     ! Read in final snapshot - specify .FALSE. to omit masses
     call read_snapshot(final_snap,.FALSE.)

     ! Apply selection criteria and find ParticleIDs - specify .TRUE. to denote this
     ! is the *original* (i.e. final snapshot) selection
     call find_selection_ParticleIDs(x0,y0,z0,r0,.TRUE.)

     ! Dump these to a master list of ParticleIDs. Specify .TRUE. to denote this
     ! is the *original* (i.e. final snapshot) selection
     call dump_selection_ParticleIDs(final_snap,.TRUE.)

  endif

  if(mode < 3) then

     ! Read in the selection - simpler coding to just assume we're using mode 2 or 3
     call read_selection_ParticleIDs(final_snap)

     ! Read in previous snapshots and cross reference to the find the indices of the selected particles. 
     ! Then compute the centre of mass of these particles, and dump along with the indices. Pass .TRUE.
     ! to read_snapshot as we now need the masses for the CoM calculation.
     write(logunit,*)
     write(logunit,*)' -finding indices of selected particles in previous snapshots...'
     call flush(logunit)
     do isnap=0,final_snap
        write(logunit,*)
        call read_snapshot(isnap,.TRUE.)
        call trace_selection_indices()
        call compute_traced_com()
        call dump_traced_com_and_indices(isnap)
     end do
     
  endif

  if(mode < 4)then 

     ! Now we want to find the ParticleIDs of all particles within r0 of the 
     ! traced centre of mass, at each epoch. 
     ! Don't need Mass this time so pass .FALSE. to read_snapshot.
     ! These are not the original selections, so pass .FALSE. to find_selection_ParticleIDs 
     ! and dump_selection_ParticleIDs
     write(logunit,*)
     write(logunit,*)' -finding all particles within r0 of the traced centre of mass'
     call flush(logunit)
     do isnap=0,final_snap
        write(csnap,'(I3.3)')isnap
        write(logunit,*)
        write(logunit,*)'  -isnap = '//trim(adjustl(csnap))
        call read_traced_com(isnap)
        call read_snapshot(isnap,.FALSE.)
        call find_selection_ParticleIDs(COM(0),COM(1),COM(2),r0,.FALSE.)  
        call dump_selection_ParticleIDs(isnap,.FALSE.)  
     end do
     
  endif

  write(logunit,*)
  write(logunit,*)' -code complete'
  call flush(logunit)

end program SelectParticles

! **********************************************************

subroutine startup()
  use ISO_FORTRAN_ENV
  use log
  use myfiles
  implicit none
  integer*4 :: nargs,iargc,ierr

  logunit = OUTPUT_UNIT

  ! Check that arguments are correct and read in parmfile
  nargs = iargc()
  if(nargs /= 7) call abort("usage: paramfile mode final_snap x0 y0 z0 r0"//char(10)//char(13)//&
       "[mode 1=all, 2=no selection, 3=no tracing]")
  call getarg(1,paramfile)
  call getarg(2,cmode)
  call getarg(3,cfinal_snap)
  call getarg(4,cx0)
  call getarg(5,cy0)
  call getarg(6,cz0)
  call getarg(7,cr0)

  read(cfinal_snap,*)final_snap
  read(cmode,*)mode
  read(cx0,*)x0
  read(cy0,*)y0
  read(cz0,*)z0
  read(cr0,*)r0

  ! Say hello
  write(logunit,*)''
  write(logunit,*)' -------------------------------------------------------------------------'
  write(logunit,*)'  SelectParticles :: Rob Crain, CAS, Swinburne'
  write(logunit,*)' -------------------------------------------------------------------------'
  write(logunit,*)'  isnap = '//trim(adjustl(cfinal_snap))		     
  write(logunit,*)'  x0    = '//trim(adjustl(cx0))		     
  write(logunit,*)'  y0    = '//trim(adjustl(cy0))		     
  write(logunit,*)'  z0    = '//trim(adjustl(cz0))		     
  write(logunit,*)'  r0    = '//trim(adjustl(cr0))		     
  write(logunit,*)' -------------------------------------------------------------------------'
  write(logunit,*)
  call flush(logunit)
  
  ! Test 
  if(mode < 1 .or. mode > 3) call abort(" unknown mode, please use 1, 2 or 3")

end subroutine startup
     
! *******************************************************************************

subroutine abort(message)
  use log
  implicit none
  character(*), intent(in) :: message
  integer*4 ::  ierr

  write (logunit,*) message
  call flush(logunit)

  stop
end subroutine abort

! --------------------------------------------------------------------
